clear all;

s = RandStream('mt19937ar','seed', 45);
RandStream.setDefaultStream(s);

dOmegaSky	= 0.1;
dOmegaEarth	= 0.05;

% Draw up a list of sky positions for the sources
skyPos		= sinusoidalMap(dOmegaSky);
nSky		= size(skyPos,1); 

% Draw up a list of Earth positions to place the detector
earthPos	= sinusoidalMap(dOmegaEarth);
nEarthPos	= size(earthPos,1);

% specify source params
% Note (Anand): even though we specify the source's Phi0, psi and iota values here, 
% we will average over these quantities when plugging into the waveform.
srcParam.Phi0   = pi*rand(1,1);
srcParam.psi    = pi*rand(1,1);
srcParam.iota   = 0;
srcParam.Delta  = pi/2 - skyPos(:,1);
srcParam.Alpha  = skyPos(:,2);
srcParam.m1     = 1.4;					%msun
srcParam.m2     = 1.4;					%msun
srcParam.mTotal = srcParam.m1 + srcParam.m2;
srcParam.eta    = srcParam.m1*srcParam.m2/srcParam.mTotal^2;
srcParam.r      = 35;					%Mpc

fixedSNR	= 10;

det.H1 = []; det.H1 = LoadDetectorInfo('H1', det.H1, 0);
[dist.H1] = calcDistAtFixedSNR(det.H1, srcParam, skyPos, fixedSNR, 'LIGO');

det.L1 = []; det.L1 = LoadDetectorInfo('L1', det.L1, 0);
[dist.L1] = calcDistAtFixedSNR(det.L1, srcParam, skyPos, fixedSNR, 'LIGO');

det.V1 = []; det.V1 = LoadDetectorInfo('V1', det.V1, 0);
[dist.V1] = calcDistAtFixedSNR(det.V1, srcParam, skyPos, fixedSNR, 'Virgo');

det.I1 = []; 

for k=1:nEarthPos
	if ~mod(k,25)
		fprintf(2, 'Cheked %6d of %6d Earth locations\n', k, size(earthPos,1))
	end
	det.I1.lat	= pi/2 - earthPos(k,1);
	det.I1.long	= earthPos(k,2);
	det.I1 = LoadDetectorInfo('I1', det.I1, 0);

	[dist.I1] = calcDistAtFixedSNR(det.I1, srcParam, skyPos, fixedSNR, 'LIGO');
	
	%-- Max distance along any given direction in the sky (3-IFO coinc)
	rmin3IFO.HLV =  min([dist.H1 dist.L1 dist.V1]');	
	rmin3IFO.HLI =  min([dist.H1 dist.L1 dist.I1]');	
	rmin3IFO.HVI =  min([dist.H1 dist.V1 dist.I1]');	
	rmin3IFO.LVI =  min([dist.L1 dist.V1 dist.I1]');

	rmin3IFO_best_HLVI = max([rmin3IFO.HLV' rmin3IFO.HLI' rmin3IFO.HVI' rmin3IFO.LVI']');

	vol3_HLV  = (1/3)*sum( (rmin3IFO.HLV).^3 )  * dOmegaSky;
	vol3_HLVI = (1/3)*sum( (rmin3IFO_best_HLVI).^3 ) * dOmegaSky;
     
	improvementMinVol3IFO(k) = vol3_HLVI ./ vol3_HLV;

end


%==== Plotting ====%
fig1 = figure(1); 
   clf; set(fig1, 'Color', 'White');

   deltaDeg	= 0.15;
   lonEarth	= 180*(earthPos(:,2))/pi;
   latEarth	= 180*(pi/2 - earthPos(:,1))/pi;
   lonGrid	= (-180:deltaDeg:180);
   latGrid	= (-90:deltaDeg:90)';
   zGrid 	= griddata(lonEarth, latEarth, improvementMinVol3IFO./max(improvementMinVol3IFO), lonGrid, latGrid);
   R		= makerefmat(lonGrid(1), latGrid(1), deltaDeg, deltaDeg);

   axesm('robinson','maplatlim',[min(latEarth) max(latEarth)],'maplonlim',[min(lonEarth) max(lonEarth)])
   gridm; 
   geoshow(zGrid, R, 'DisplayType', 'texture'); 
   load coast;
   geoshow(lat, long, 'Color', [0.15 0.15 0.15]);

   mlabel('MlabelLocation',[-180 -135 -90 -45 45 90 135 180], 'MLabelParallel', 'equator', 'FontSize', 10);
   plabel('PlabelLocation',[-90 -60 -30 30 60 90], 'FontSize', 10);
   colormap jet;
   cb = colorbar;
   xlabel(cb, 'Location Merit');
   caxis([0.45 1]);

   % Plot the location of all the detectors
   plotm(det.V1.lat*180/pi, det.V1.long*180/pi, 'ro', 'MarkerSize', 6, 'MarkerFaceColor', 'w'); 
   plotm(det.L1.lat*180/pi, det.L1.long*180/pi, 'ro', 'MarkerSize', 6, 'MarkerFaceColor', 'w'); 
   plotm(det.H1.lat*180/pi, det.H1.long*180/pi, 'ro', 'MarkerSize', 6, 'MarkerFaceColor', 'w'); 

   tightmap;
   print -dpng 'LocationMerit_intialLIGO.png'




